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Abstract 

We present a general scheme for treating the integrable singular terms within exact exchange 
(EXX) Kohn-Sham or Hartree-Fock (HF) methods for periodic solids. We show that the singularity 
corrections for treating these divergencies depend only on the total number and the positions of 
k points and on the lattice vectors, in particular the unit cell volume, but not on the particular 
positions of atoms within the unit cell. The method proposed here to treat the singularities 
constitutes a stable, simple to implement, and general scheme that can be applied to systems with 
arbitrary lattice parameters within either the EXX Kohn-Sham or the HF formalism. We apply the 
singularity correction to a typical symmetric structure, diamond, and to a more general structure, 
trans-p olyacetylene. We consider the effect of the singularity corrections on volume optimisations 
and k point convergence. While the singularity corrections clearly depends on the total number of 
k points, it exhibits a remarkably small dependence upon the choice of the specific arrangement of 
the k points. 
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I. INTRODUCTION 



In recent years exact exchange (EXX) Kohn-Sham (KS) methods for sohds became in- 
creasingly popular-i^ii^i^'^i^ as alternative to conventional KS procedures based on the local 
density approximation^!^ (LDA) or generalised gradient approximations (GGAs).- EXX-KS 
methods treat both the exchange energy as well as the local KS exchange potential, not to 
be confused with the non-local Hartree-Fock exchange potential, exactly. This means they 
constitute a systematic improvement over LDA and GGA methods in the sense that only the 
correlation energy and potential, i.e., contributions of higher order in the electron-electron 
interaction, need to be approximated whereas the terms of first order in the square of the 
electron charge, i.e., the Coulomb and exchange energy and potential, are treated exactly.— 
Because exchange and Coulomb potential and energy are treated exactly unphysical self- 
interactions contained in the Coulomb energy and potential are completely cancelled by the 
exchange energy and potential. EXX methods therefore are free of Coulomb self-interactions. 
As a result EXX band structures and in particular band gaps are strongly improved com- 
pared to those from LDA or GGA methods. Indeed for medium gap semiconductors EXX 
methods yield band gapa^ which are very close to the experimental ones^i despite the fact 
that the correlation potential needs to be neglected or approximated by conventional LDA 
or GGA functionals, and despite the fact that the KS band gap does not account for the 
derivative discontinuity— of the band gap at integer electron numbers. 

A second first-principles approach besides the family of density-functional methods is 
the Hartree-Fock (HF) method.— Recently, there has been an increasing interest in HF 
methods for solids as basis for higher level approaches like, e.g., M0ller-Plesset,— coupled 
cluster— or multireference configuration interaction^^ methods. 

Both in the EXX and in the HF formalism the exchange energy contains divergent terms.— 
In the limit of an infinite system, i.e., the limit of an infinite number of k points, the diver- 
gencies are integrable. In this limit the exchange energy is therefore well-defined. Moreover, 
corresponding divergencies also occur in the matrix elements of the non-local exchange oper- 
ator which is required in the HF self-consistency process and can be used in the construction 
of the local KS exchange potential.-*^ Also here the divergencies are integrable in the limit 
of an infinite number of k points. The question arises how to treat these divergencies in 
practical calculations which necessarily take into account only a finite number of k points. 
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Indeed, in order to keep the computational effort as low as possible, it is preferable to keep 
the number of k points as low as possible. This, however, is possible only if an adequate 
treatment of the singularities is available. Moreover, such a treatment of the singulari- 
ties should be computationally efficient and ideally its implementation should not require 
much programming effort. Gygi and Baldereschi^ presented such a method for the case of 
zincblende (fee) structures. Wenzien et al.— further generalised the method to simple cubic, 
bcc, hexagonal, and orthorhombic structures. For other crystal structures such simple and 
straightforward method, to our knowledge, is still lacking and alternative approaches^>22ii^ 



are more involved. In Refs. 
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]22[ |. e.g., a general treatment of the singularities is presented. 
This method, however, is somewhat laborious because it requires a quadrature over recip- 
rocal lattice vectors at each k points. This quadrature formally has to run over an infinite 
number of reciprocal lattice vectors which in practice needs to be approximated by a finite 
summation. Thus there is demand for a simple, efficient treatment of the singularities, that 
is applicable to arbitrary crystal structures. 

In this article we present a simple, efficient, and general treatment of the singularities 
in Hartree-Fock and exact-exchange Kohn-Sham methods for periodic systems, which ex- 
tends the approach of Gygi and Baldereschi^° to systems with arbitrary lattice parameters. 
The derivation of this treatment of the singularities is accompanied by an analysis of the 
singularities and demonstrates the simplicity of the suggested method for handling these sin- 
gularities. In order to demonstrate the applicability of the approach we present results for 
the symmetric diamond (fee) structure (2 carbon atoms) as well as for trans-polyacetylene 
(4 carbons and 4 hydrogens in its crystalline unit cell). Polyacetylene has monoclinic sym- 
metry P2i/a, i.e., non-orthogonal lattices,— and constitutes a simple example of an organic 
polymer. 

The article is organized as follows. In Section [III the general treatment of the divergent 
terms in EXX and HF methods is derived and discussed. Section IIIII unwraps the results 
for diamond and trans-polyacetylene. Section [IV] concludes. 

II. TOTAL ENERGY WITHIN THE EXX FORMALISM 

We start by considering the expressions for the total electronic energy Eq within the 
Kohn-Sham (KS) and the Hartree-Fock (HF) schemes. Within the KS formalism the total 
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ground state energy is decomposed into the non- interacting kinetic energy Ts, the Coulomb 
energy U, the exchange energy E^, the correlation energy Ec, and the interaction energy 
with the external potential, V^xt'^ 

Eo=Ts + U + E, + E, + jdr K.t(r)po(r). (1) 

The non-interacting kinetic energy T<j is evaluated exactly via the KS orbitals. The contribu- 
tions U and j drVf-xti'^) Pq{t^) ^Iso can be calculated exactly for a given electron density and 
thus also for the ground state electron density po- The correlation energy E^ in almost all 
KS schemes is evaluated approximately within the LDA~ or the GGA.- The exchange energy 
Ex can either be evaluated via the LDA or the GGA within a conventional KS scheme, or 
exactly within the EXX-KS scheme.— 

The HF total energy, on the other hand, is decomposed into 

Eo=T + U + Ex + jdv Vext{r)pHF{r). (2) 

Within HF schemes all contributions of the energy are usually treated exactly: the kinetic 
energy T and the exchange energy Ex via the HF orbitals, and U and J dTVext{'^)pHFij) via 
the HF electron density Phf- 

The exact-exchange energy Ex per unit cell for a crystalline solid, either for the KS or 
HF formalisms, is given by: 

Dk wq |r r I 

where both summations run through all occupied single-particle wave functions, i.e., orbitals 
0^,k and for each k point in the Brillouin zone (BZ). All orbitals are assumed to be 
normalised with respect to the crystal volume VL = N^V where V designates the volume 
of the unit cell, and Nk denotes the number of k points. We implicitely treat the spin 
via appropriate prefactors in summations and consider for simplicity non-spin polarised 
calculations. The Coulomb interaction term, j^r^ i'^ Eq. ([3]), has to be treated taking into 
account periodic boundary conditions. Note that, despite the fact that the expression for the 
exchange energy in terms of one-particle functions is identical in the KS and HF case, the KS 
and HF exchange energies remain different because their respective one-particle functions 
are constructed using two different scheme: KS uses a local exchange operator while HF 
uses a non-local exchange operator. 



After expressing the product of one-particle functions as 



^U(r)0.k(r) = ^Y1 >^-c,,.k(G)e^(«+'^-'^)-, (4) 

G 



with 



Y^^MG) = Idv e-(«+i'-'i)- 0l,^(r) 0,k(r) (5) 
Jn 

one obtains the following expression for the exchange energy per unit cell 

A-K ^u^q,t,k(G) ^toq,i,k(G) 

^^ = -N^l^l^l^ |G + k-q|2 ' 

if the following relation is taken into account 



g-iG-rgiG-r' 

w 



I — 1:^ — :7i — = J7m ^gg'' C^) 



which holds due to translational symmetry. 

Expression ([6]) contains singular terms, namely those with G = 0, k = q, and v = w. 
Note that when f 7^ w no singularities occur for any value of G and k. This is due to the 
relation 

ywk,vk{0) = Swv (8) 

which holds because Eq. ([5]) that defines Ft„q,^k(0) in the case where G = and k = q just 
turns into the orthonormality condition for the one-particle functions. Thus, contributions 
with G = 0, k = q, and v ^ w vanish because the plane wave representations of the products 
'i^ik(^)*^fk(i^) with V ^ w do not contain contributions from a plane wave with G = 0. [This 
means that for v ^ w no singularities are present in Eq. (I6l). Therefore, stricktly speaking, 
Eq. ([H]) needs to be modified in a way that for v ^ w singular terms are no longer present]. 

Due to the presence of the singularities described above the exchange energy is well- 
defined only in the limit of an infinite number of unit cells, i.e., for — > 00. In this case, 
the summations over k and q turn into integrals over the BZ and Eq. ([6]) for the exchange 
energy assumes the form 

47r VI? f ^ f ^«^q,Dk(G) V'«;q,i,k(G) 

The singularities in integral ([9]) are integrable. Therefore, the exchange energy is now well- 
defined. 
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Adopting an idea of Gygi and Baldereschi^S we now manipulate the contributions on the 
right hand side of Eq. iQ with G = and v = w, i.e., those contributions which contain 
the integrable singularities, by adding and subtracting a function /(q) which shall obey the 
three following conditions: (i) /(q) is periodic within the reciprocal lattice, (ii) /(q) diverges 
as for q ^ and is smooth elsewhere, and (iii) /(q) = /(— q)- This leads to 



' ~Nkn (27r)6 
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and 
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The function Fy_ in Eq. (fTT!) is given by 



47r 



^k^^E/(k-^)- 



(10) 



(11) 



(12) 



(13) 



qy^k 



In Eq. f ITU]) A^t, designates the number of valence bands, which comes from the summation 
over the valence bands in the two terms containing the function /. We also used condition 
(ii) for the function / and Eq. ([8]) that require Y*^^^J^Q) 'K„q,i,k(0)/|k - qp - /(k - q) for 
any given k to be a smooth function of q that equals zero at q = k. Therefore, the first 
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integral over q and k after the first equality sign of Eq. (fTOi) can be evaluated by summations 
over the finite grid of k points omitting the terms with q = k. Due to the periodicity and 
inversion symmetry of /(q) [condition (i) and (iii) for /] the integrals of /(k-q) over the BZ 
can be replaced by BZ integrals of /(q). Furthermore, for the case of a uniform grid of k 
points the function F of Eq. fllip simplifies to 

^ = fE/(q)- (14) 

q^O 

The evaluation of the exchange energy can now be done according to 

4^ ^ ^ Y-^4,.k(G)K.q,.k(G) 



iVfcfi^ r |G + k-q|2 

+N,[F-F]. (15) 

This implies that for the evaluation of the exchange energy the singular terms in the original 
expression ([6]) can first simply be omitted and then be taken into account by N^[F — F], 
i.e., by adding Ny[F — F] to the exchange energy obtained if the singular terms are simply 
omitted. The correction is calculated only once before the self-consistency procedure. In 
fact, the correction N^[F — F] depends only on the unit cell lattice vectors, and thus in 
particular on the unit cell volume V , and on the number A^^^ and the positions of the k 
points. It does not depend on the number, type, or positions of the atoms within the unit 
cell, and does not depend on the one-particle wave functions. This has obvious advantages 
for atomic relaxations at fixed unit cell volumes and fixed lattices. 

The whole scheme, of course, hinges on the availability of a suitable function /(q). For 
fee systems such a function was given by Gygi and Baldereschi.-^ For sc, bcc, hexagonal, and 
orthorhombic systems Wenzien^i presented such functions. Here we suggest the following 
function / for arbitrary crystal structures: 



-1 



+ 2 5]][bi sin(aj- ■ q)] ■ [b^+i sin(aj+i ■ q)] I . (16) 
i=i J 

The bj (with = bi for a compact formulation accounting for cyclic permutations) are the 
reciprocal lattice vectors, and the a^ (with = ai) are the corresponding lattice vectors 



spanning the unit cell. The coefficient l/(27r)^ arises from the factor 27r contained in the 
Taylor expansion of the trigonometric functions, because a^-q implicitely contains aj hj = 27r 
if q is expressed as q = with qj describing the components of q with respect to 

reciprocal lattice vectors. The function flTB]) by construction has the required periodicity 
of the reciprocal unit cell. Expansion into a Taylor series with respect to the cartesian 
components Qx, Qy, and g^, or equivalently with respect to qi, q2, and q^, the components of 
q referring to the reciprocal lattice, furthermore shows that it diverges as for q ^ 0. 
Therefore, /(q) satisfies requirements (i)-(iii) for any type of (linearly independent) lattice 
parameters: ai, sl2, and as. 

The integration over the BZ of the function (1161) required for obtaining the correction 
F in Eq. f[T^ can easily be carried out numerically on an adaptive grid using an iterative 
algorithm. To this end, we place the reciprocal lattice centered symmetrically around q = 0. 
In the ffist iteration we generate a regular (2A^ + 1) x [2N + 1) x (2A^ + l)-grid with the 
number N being a multiple of 3, typically = 60.— The grid points shall be labeled q^m„ 
with -N < i < N, -N < m < N, and -A^ < n < N, with the point qooo located at 
the origin of the reciprocal lattice. We then divide the unit cell into an inner part given by 
a cell in reciprocal space which again is centered symmetrically around q = 0, and which 
is defined by lattice vectors being one third of the original reciprocal unit cell vectors. In 
the first iteration numerical integration is carried out only in the outer region. Then in the 
inner region, where the singularity is located, the number of points is tripled and a second 
iteration proceeds as the ffist one, working now on the outer part of the inner region of the 
first iteration. By moving on in this fashion the algorithm triples the mesh size around the 
singularity in each iteration step. Therefore, the integration result is more accurate than 
with any regular mesh. We observe that less than 10 steps, depending on the lattice vectors 
considered, are sufficient to get the integral converged. The implementation of the described 
numerical integration is straightforward leading to about 200 lines of Fortran instructions. 
The computational time for carrying out the integration is negligible. 

Having considered in detail the treatment of the singularities in the exchange energy 
we now briefly present the corresponding treatment of singularities in the evaluation of 
the matrix elements of the non-local exchange potential, which is required in the HF self- 
consistency process, or can be used in the construction of the local KS exchange potential^i^ 
during the self-consistency process of a KS calculation. 
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The matrix elements of the non-local exchange potential, v^^(k,^, u), are given by 

NLn \ /■ , , ' XMk(r')0»q(r')0i,q(r)Xz.k(r) 

(k, /U, z/) = - 2^ / drdr , _^,| (17) 

with Xfj-k and Xuk denoting the basis functions for the representation of the one-particle 
functions (pwci- The basis functions x^k ^.re products of a periodic part and a Bloch factor 
g-jk r rjj^g most common choice for the basis functions x^k are plane waves e~*('^+'') '". Like 
in the treatment of the exchange energy, we now express the products 0|„q(r) Xi^k(r) as 



( 
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Expression fll9p contains singular terms, again those with G = and k = q. In the limit of 
an infinite number of unit cells the summation over q again turns into an integral, namely 
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with an integrable singularity. We can now treat the singular terms in the right hand side 
of (120|) exactly analogously to the singular terms occuring in the exchange energy: 
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Thus, the matrix elements f^'(k, /i, i/) of the non-local exchange potential can be eval- 
uated by first omitting the singular terms in Eq. f[T9l) and by then adding the following 
correction term 



E^-Vk(0)Wk(0) 



(22) 



The required sums Fk and the integral F have to be calculated only once at the beginning 
of the self-consistency procedure and then are multiplied by ^ujk,/ik(0) ^ti)k,!/k(0) in each 
HF self-consistency cycle, because the y^k,^tk(0) ^«ik,!/k(0) change during the self-consistency 
cycle. In case of a uniform grid of k points the reduce to F given in Eq. (1141) . 

The most widely used basis set for the one-particle functions of periodic systems are plane 
waves. If basis sets of plane waves are employed then the one-particle functions 0«,q(r) are 
given by 

(23) 



-i(G+k)-r 



The elements f[T9|l of the non-local exchange potential turn into 
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and the correction term fl22D turns into 
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(25) 



Again the matrix elements v^^(k, G, G ) can be calculated by first omitting the singular 
terms in Eq. fl24l) and by then adding the correction term (l25!l . 

So far we have considered only systems with bands that are either fully occupied or fully 
unoccupied., i.e. we have considered isolating systems at zero temperature. In the Appendix 
we sketch how the formulas of this Section change for systems with partially filled bands 
and how the singularities can be treated in this case. 



III. RESULTS FOR DIAMOND AND Ti? A iVS-POLYACETYLENE 

The applicability of the presented approach for treating the divergencies is now demon- 
strated by applying it to two cases, diamond and trans-polyacetylene, using different approx- 
imations for the exchange-correlation functionals: the Slater-Dirac (exchange-only LDA) re- 
ferred to as Dirac exchange in the following, the complete exchange-correlation LDA in the 
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parametrisation of Vosko, Wilk, and Nusair (VWN),- the combination of Dirac exchange 
plus Perdew86^ (P86) correlation (P86 being VWN correlation plus a gradient correction), 
the EXX (exact-exchange only), and finally the combination of EXX with P86 correlation. 

The pseudopotentials were generated using the pseudopotential generation code of 
Engel,-- which is based on the TrouUier- Mart ins norm-conserving scheme.— In all cases 
the pseudopotentials were generated using consistently the same functionals for exchange 
and correlation as for the plane wave calculations. Relativistic effect are not included. The 
pseudopotentials of C are all constructed using a cutoff radius of 1.3 a.u. for both s and p 
levels. We constructed for the calculations of trans-polyacetylene chains hydrogen pseudopo- 
tentials with a cutoff radius of 0.9 a.u.. For diamond, the energy cutoff of the plane-wave 
basis is 60 Ry for the one-particle functions, and 20 Ry for the exchange potential and the 
response function^. For trans-polyacetylene,— we reduced the energy cutoffs to 32 Ry for 
the one-particle functions and 12 Ry for the exchange potential and the response function, 
because we are interested to determine the effect of the singularity function in terms of 
several possible k points meshes, and study meshes with a large number of k points. 

The lattice constants for diamond were varied from 3.1 to 4.1 A. For comparison, the 
experimental lattice constant of diamond is 3.5668A.'ii Figure [1] shows the variation of the 
singularity correction of the EXX exchange energy as a function of the k point mesh for 
diamond. Fig. [2] shows the singularity correction for diamond for a fixed number of k 
points but for different volumes. Choosing 5 x 5 x 5 k points ensures convergence of 
the total energy within 0.2 eV, while 8 x 8 x 8 k points ensure total energy convergence 
within 0.05 eV. From Fig. [1] we notice that the singularity correction Ny{F — F) and the 
exchange energy without the singularity correction vary oppositely with increasing number 
of k points. The complete exchange energy including the singularity correction turns out to 
be quite stable with the number of k points. Fig. [1] also shows that for small and medium 
numbers of k points a more symmetric mesh division reduces the deviation of the complete 
exchange energy from its converged value at high numbers of k points (compare dashed with 
continuous EXX line). Figure [2] also shows an aspect important for volume optimisations: 
the singularity correction changes dramatically with the volume and thus strongly modifies 
the position of the energy minimum as well as the bulk modulus. The energy minimum is 
reduced because the correction function is monotonically increasing with the volume. The 
bulk modulus is also modified because the variation of the correction is obviously not linear. 
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The bulk modulus of diamond without the singularity correction is much smaller (25% 
smaller) than with it. Therefore, accurate integration of the singularity in the exchange 
energy is essential for evaluating bulk properties within EXX or HF methods. 

Fig. [3] shows volume optimisation results for diamond using the various combinations of 
exchange and correlation functionals. We observe that a removal of Coulomb self-interactions 
(see EXX versus Dirac exchange) induces a significant reduction of the total energy (~2 eV). 
It also leads to a reduction of the lattice constant minimum [compare vertical lines in Fig. [3] . 
The values of both exchange-only energy curves, i.e., of the EXX and Dirac-Slater curves, are 
much higher than the curves that contain a correlation potential, i.e., the LDA, Dirac+P86, 
and EXX+P86 curves, which reflects that correlation affects the total energy. The reduction 
of the total energy from EXX to EXX-P86 is of the same order as the reduction from Dirac- 
exchange only to the Dirac exchange plus P86 correlation. However, the lattice constant 
minimum of the EXX-P86 is shifted to a much lower value than any other of the combinations 
of functionals. The reason for this maybe the reintroduction of self-interations, through the 
P86 correlation function. In any case, the poor performance of the combination EXX-P86 is 
not surprising because the P86 correlation is not meant to be used with the EXX, but rather 
with the LDA or GGA exchange in order to exploit error cancelations between exchange 
and correlation. Therefore, development of correlation functionals that do not depend on 
such error cancelations and thus are well-suited for combination with the EXX is highly 
desirable. 

The EXX energy optimised lattice equals 3.555 A [see Fig. [3]. It becomes natural now to 
evaluate the band structure at the EXX energy minimum. Figure H] shows the band structure 
at this EXX energy minimum. The indirect EXX bandgap is 4.838 A, a value comparable 
to previous published data,^ and much closer to the experimental band gap of 5.50 eV— 
than the LDA value of 3.90 eV.^- The experimental lattice constant of diamoncUi (3.5668 
A) is slightly larger (+0.011 A) than the EXX energy optimised lattice. Going from the 
EXX lattice minimum to the experimental value, i.e., addition of 0.011 A to the EXX lattice 
constant, leads to minute reduction of the band gap. However, we want to emphasize that the 
singularity correction to the exchange energy is obviously essential for determining the right 
correspondence between the EXX energy lattice minimum [Fig. [3] and the band structure or 
its band gap [Fig.H]. For instance, without the singularity correction the EXX lattice energy 
minimum would be incorrectly overestimated (3.730 A instead of 3.555 A as illustrated in 
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Fig. [3]), and the EXX band gap correspondingly would be largely underestimated (because 
the band gap variation is roughly inversely proportional to the lattice constant). In summary, 
we find that for diamond the EXX band gap is somewhat smaller than in the experiments 
(0.66 eV lower than experiment), but EXX improves significantly the LDA value (~1.6 eV 
lower than experiment). As mentioned in the introduction, the EXX calculation does not 
account for the derivative discontinuit}*^'^^ of the band gap at integer electron numbers and, 
of course, also not for the correlation potential. 

We now consider a more general crystal structure: trans-polyacetylene. The unit cell 
is constituted of 4 carbons and 4 hydrogens. More data on the band structure of trans- 
polyacetylene can be found elsewhere.— Trans-polyacetylene constitutes a monoclinic lat- 
tice structure (group Pixjd) with the following lattice parameters expressed in cartesian 
coordinates (and in A): 

ai = (4.24,0.00,0.00) 

a2 = (-0.0642644,2.454158,0.00) (26) 
aa = (0.00,0.00,7.32) 

The angle between ai and a2 is 91.46°. The angle between any two dimerised chains is 
55°. The coordinates of the structurally optimised hydrogen atoms come from Hartree- 
Fock calculations and the lattice parameters and C-C bond distances and angles come from 
experimental values.— This structure constitutes a general and realistic case to test our 
singularity correction for several k points meshes. A graph equivalent to Fig. [1] is displayed 
in Fig. [5] for this molecular crystal. 

Figure [5] shows for trans-polyacetylene a similar trend as shown in Fig. [1] for diamond. 
That is, the singularity correction and the exchange energy excluding the singularity vary 
in opposite matter, for any chosen k points division. The complete exchange energy is a 
rather monotonic function of the total number of k points in the unit cell. 

The results, both for diamond and trans-polyacetylene, show that the approach presented 
here to treat the integrable singularities in the KS and HF methods constitutes a stable and 
general scheme. 
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IV. CONCLUDING REMARKS 



We have presented a general scheme for treating the integrable singularities of the ex- 
change energy within the EXX or the HF formahsms. We have shown that the divergent 
terms in the exchange energy depend only on the number and positions of k points and on 
the unit cell vectors and thus on the unit cell volume, but not on the single particle wave 
functions or on the particular atomic positions within the unit cell. A similar correction 
procedure is proposed for matrix elements of the non-local exchange operator which occurs 
in the Hartree-Fock methods and can be used to construct the exact local Kohn-Sham ex- 
change potential. We applied the singularity correction to a typical symmetric structure, 
diamond, and to a more general structure, ^rans-polyacetylene, and discussed the effect of 
the singularity function on volume optimisation and k points convergence. The singularity 
function depends strongly on the total number of k points and more weakly on the choice of 
the specific division of the k points mesh. The complete, i.e., singularity corrected exchange 
energy, converges well with the number of k points. The method proposed here constitutes 
a stable, simple to implement, and general scheme that can be applied to systems with any 
lattice parameters within either the EXX Kohn-Sham or the Hartree-Fock formalism. 
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VI. APPENDIX 



In this Appendix we briefly consider the treatment of singularities for systems with par- 
tially filled bands. In this case expression for the exchange energy turns into 
1 o^-occ ^ /■ , 0L(r)0^q(r)0i,q(r>,i,(r) 

J^k Jo Jo r — r 

with the occupation factor r/^k. The occupation factor shall not include the spin multiplicity, 
i.e., < ?7t,k < 1. Analogously to Eq. ([6]) the exchange energy can be expressed by 

= - L L ^^'^ L |G + k-qP ' (2^) 

The singular terms in expression ( l28l) . namely those with G = 0, k = q, and v = w, can be 
treated in analogy to Eq. (ITUl) according to 
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with Fk = f Eq^k/(q - k) and F = f^^dq /(k - q). In Eq. m we have used 

that the integral j^^dkrjl^, does not contain any singularities and therefore can be evaluated 
via summation over the k points. The energy again can be evaluated by first omitting the 
singular terms in Eq. ( 128|) and by then adding the correction term 



.2 
/i;k 



F. 



(30) 



For the particular case of a uniform grid of k points the functions Fk all equal the function 
F of Eq. (fT2|) and the correction term turns into 



[F-F] 



(31) 



In a similar way as the treatment of the singularities in the exchange energy was gener- 
alized for the case of partially occupied bands also the treatment of the singularities in the 
exchange potential can be generalized to the case of partially occupied bands. The matrix 
elements of the non-local exchange potential, f^'^(k, /x, z/), then are given by 



4:71 ^«^q,/ik(G) ^u.q,i/k 



|G + k-q| 



(32) 



uiq G 

Expression (!32|) contains singular terms, i.e., those with G = and k = q. In the limit of 
an infinite number of unit cells the summation over q again turns into an integral, namely 

47r n 

' BZ 



LJ^V..!^ |G + k-qP 



(33) 



with an integrable singularity. We can now treat the singular terms in the right hand side 
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of (133|) exactly analogously to the singular terms occuring in the exchange energy: 
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I K — 111" 

qT^k 

Thus, the matrix elements f^'(k, /i, i/) of the non-local exchange potential can be eval- 
uated by first omitting the singular terms in Eq. fl521) and by then adding the following 
correction term 



E^-kl^\Mk(0)n,k,.k(0) 



F . 



(35) 



The required sums Fk and the integral F have to be calculated only once at the beginning 
of the self-consistency procedure and then multiplied by Xlw, ^^^k ^u^k,/ik(0) ^uik,i^k(0) in each 
HF self-consistency cycle, because the ^ujk,/ik(0) ^i«k,i/k(0) and the occupation numbers rj^k 
change during the self-consistency cycle. In case of a uniform grid of k points the Fk reduce 
to F given in Eq. f fT4l) . as for fully occupied bands (since the singularity function /(q) does 
not depend on the occupation numbers r^^k)- 

If basis sets of plane waves are employed then the elements flH21) of the non-local exchange 
potential turn into 



<nk,G,G 



An Cy^diG — G + G )C*{G ) 



n 



toq q" 

and the correction term fl35|) turns into 



VwVi C^k(G)C*k(G 



|G - G +k-q| 



(36) 



(37) 



Again the matrix elements f^^(k, G, G ) can be calculated by first omitting the singular 
terms in Eq. (l36l) and by then adding the correction term (1371) . 
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# k points (log scale) 



FIG. 1: Singularity correction of the exchange energy as a function of the number of k points. 
Circles represent data for the EXX energy without taking into account the singular terms. Triangles 
are the singularity correction. Squares constitute the full EXX energies, including the singularity 
correction. The exchange energy excluding singular terms and the singularity correction vary 
oppositely with the number of k points and their sum leads to a relatively constant and quite fast 
converging EXX energy. The lines guide the eyes for non-symmetric (dash) and symmetric (solid) 
k points meshes. The number of k points along axes of the reciprocal lattice are indicated by the 
numbers in parenthesis. 



20 




7 8 



10 11 12 13 14 15 16 17 18 

° 3 

Volume [A ] 



FIG. 2: Singularity correction of the exchange energy of diamond for fixed k point mesh (5x5 
X 5) as a function of the volume. [See text for details]. 
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FIG. 3: Volume optimisation for diamond. We used the equation of states of Teter.-^i The total 
energy as a function of volume is depicted, using different exchange-correlation functionals: the 
Dirac exchange (exchange only LDA), LDA, Dirac exchange plus P86 correlation, EXX with or 
without P86 correlation, and EXX without singularity correction. The singularity correction, as 
depicted in Fig. [21 leads to a significant shift of the energy- volume minimum for EXX, corresponding 
to a reduction of the lattice constant of -0.175 A. (Compare curves with stars and open-diamonds). 
For a discussion of correlation effects with EXX-I-P86 and Dirac-|-P86, see text. The experimental 
lattice constant of diamond is 3.5668 A (Volume = 11.3443 A^). 
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FIG. 4: Band structure of diamond evaluated at the EXX optimised lattice constant of 3.555 
A[See Fig. [3]. The band gap of diamond is indirect, towards the T~X direction of the BZ. The 
EXX direct transition at F equals 6.253 eV. The EXX band gap equals 4.738 eV and is located at 
72% of the X point away from F. The experimental lattice constant and band gap of diamond are 
respectively 3.5668 A and 5.50 eVjii 
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FIG. 5: Singularity correction of EXX exchange energy of irans-polyacetylene as a function of the 
number of k points. See caption of Fig. [T]for symbols description. The number of k points along 
axes of the reciprocal lattice are indicated by the numbers in parenthesis. The second entry refers 
to the number of k points along the irans-polyacetylene chain. 
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